function model = randInitial(nos, k, r, param)

[u,s,v] = svd(nos,0);
[p,m] = size(nos);
if (~isfield(param,'U0'))
    U = u(:,1:r) * sqrt(s(1:r,1:r));
    V = (sqrt(s(1:r,1:r)) * v(:,1:r)')';
else
    U = param.U0;
    V = param.V0;
end

Sigma = zeros(1,1,k);
for num = 1:k
    Sigma(:,:,num) = rand;
end

R = initialR(nos(:)', k);
nk = sum(R,1);
w = nk/(p*m);

model.U = U;
model.V = V;
model.Sigma = Sigma;
model.weight = w;
model.k = k;

model = Maximization(nos, R, model,0);

function R = initialR(X, init)
[d,n] = size(X);
k = init;
idx = randsample(n,k);
m = X(:,idx);
[~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
[u,~,label] = unique(label);
while k ~= length(u)
    idx = randsample(n,k);
    m = X(:,idx);
    [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
    [u,~,label] = unique(label);
end
R = full(sparse(1:n,label,1,n,k,n));
